I.  Introduction

II.  Four-Scale model

Trees are modelled as discrete objects; for conifers 4-Scale uses a cone plus a cylinder shape [27] (see Fig. 1) to represent the crown. The crown's surface is treated as a complex medium where shoots (or leaves for deciduous canopies) can cast shadows on each other so that even on the sunlit side of a crown, the viewer can see shaded foliage. The patchiness of the stand tree distribution is simulated with the Neyman type A distribution [25]. The Neyman distribution uses the probability of having groups of trees randomly distributed in a forest. Both the group size and the location of the group follow a Poisson distribution, hence the notation double Poisson distribution is also used to describe it. The gap fraction is computed using the Neyman distribution with an averaged gap probability within the tree crowns and a negative binomial function. These functions allow the sun's beam and view line penetration into multiple trees to compute the overlapping of the crowns. The hotspot is modelled with two kernels, one uses the gap size distribution within the crowns for the canopy hotspot and the other one uses the gap size distribution between the tree crowns for the ground hotspot. The size of the leaves or shoots, the density of tree stands, and the crowns size are important factors in determining the shape of the hotspot. 

The model parameters can be separated in three categories; 
1) sites parameters (domain size, leaf area index, tree density, tree grouping index with quadrat size and solar zenith angle), 

2) tree architecture parameters (crown radius and height, apex angle, needle-to-shoot ratio, foliage clumping index and tree foliage typical size), and

3) reflectivities of the foliage and the ground. These later parameters depend on the sensor and the wavelength used and thus can be seen as more adjustable than the architectural parameters. 

The model finally computes four surface proportions: foliage illuminated (P_T), foliage shaded (Z_T), and ground illuminated (P_G) and shaded (Z_G). Each proportion is multiplied by its reflectivity factor that depends on the wavelength used: 

A.  Branch architecture

The probability of j branch overlapping in the direction of the beam over a length `s(q) can be expressed using the Poisson distribution as [6]:

In these branch architecture simulations, the length of the branches must not exceed the height of the crown as a_b gets close to 90^°. For large values of a_b , it is preferable to use a fix value for R_b.

B  Spheroid crown shape

The mutual shadowing from one crown to another used in the 4-Scale model for simulating the proportion of sunlit crown with height was approximated by dividing the crown into two components: the cone and the cylinder and with the penetration through multiple crowns with Q_1tot and Q_2tot. For the spheroid shape, the mutual shadowing effect is only considered with Q_1tot and Q_2tot with a minimal effect on the final proportions since the division of the crown into two components (sunlit and shaded) for the conifers serves only as a weighting factor for Q_1tot and Q_2tot.

C  Overlapping of tree crowns with the repulsion effect

D.  Bi-layer forest

Assuming that the lower canopy is much shorter than the upper one, the upper canopy is then independent of the lower one, i.e., the shadow cast by the lower canopy doesn't affect the view and illumination of the upper canopy. The new process can be separated in three steps:

Step 1: the original 4-Scale model is run with the parameters from the upper canopy, a subscript 1 is added to the notation for the upper canopy, and subscript 2 for the lower canopy. The only parameter change is the effective height, which is used in the calculation of the ground hotspot kernel, that must be reduced to take into account the lower canopy. This is done by reducing the value of H_a1 (the height of the part without foliage of the upper canopy) by the total height of the lower canopy (H_a2+H_b2 + H_c2). The change in H_a1 affects only the in-between crowns hotspot calculation. Because the trees are assumed to be at equal height, the overshadowing of tree crowns does not depend on the crown base height.

Step 2: the original 4-Scale model is run with the parameters of the lower canopy. It must be noted that for the hazelnut canopy, the overlapping of the crowns is allowed at nadir to better represents the real canopy.

Step 3: with the four components of step 1 and 2, plus the gap fraction from both canopies, the six portions of the double-canopy structure can be computed.

The two components of the upper canopy (P_T1 and Z_T1) are taken directly from step 1. The probability of seeing illuminated foliage on the lower canopy if the upper canopy is absent was computed in step 2 (P_T2), but the solar beam that reaches the lower canopy and is reflected to the viewer has the probability of P_G1, thus the probability of seeing illuminated foliage in the lower canopy in the presence of an upper canopy is simply

III Comparison between model and airborne 
POLDER data

A. Sites description

1)  Old black spruce

2) Young jack pine

3)  Old jack pine

4)  Old aspen

B. POLDER data

1) The POLDER instrument

2) Data acquisition

The POLDER concept also allows the derivation of directional measurements on each grid point of extended areas surrounding the tower site. This has been applied to the datasets of Old Black Spruce and Young Jack Pine, acquired on July 21, for an area of 5 ×5 km² with a spatial resolution degraded to 100 ×100 m². For a sharp investigation of forest directional reflectance, the ``typical" BRDF of POLDER, as shown in Fig. 7c and 7d, does not use the optimal angular resolution that can be achieved by POLDER. For both datasets of Old Black Spruce and Young Jack Pine, we suppose that the adjacent pixels of the supersites are identical to those of the supersites within an area of 900 ×900 m². This hypothesis is judged realistic because it is known from available cover maps that all the pixels within the window of 81 pixels belong to the same species. All these pixels have a directional sampling slightly different from each other because they are observed under different angular conditions. Therefore, collecting BRDF measurements on points within 900 meters from the tower sites and superimposing them, we compose a new BRDF with a high directional resolution. As about 40 reflectance measurements per pixel are obtained, a complete BRDF with 3240 points (81 pixels × 40 reflectance) is thus produced from both sites (Figs. 7a and 7b). Although the 81 pixels in each site seem to have quite a large tree height variation (7.5 m-12.5 m for Old Black Spruce, 2.5 m-7.5 m for Young Jack Pine)., the dominant trees are more similar. The height variation, not included in the model, can underestimate the amount of shaded foliage in the smaller trees due to the shadow cast by taller trees.

The atmospheric correction algorithm 6S [34] is applied to each measurement. The aerosol optical depth of the full atmosphere and the fraction below the aircraft, obtained from the BORIS database, were interpolated at 550 nm and used as inputs to the algorithm. The following values were respectively deduced for the total and below-aircraft optical depths Old Aspen (0.115, 0.075), Young Jack Pine (0.115, 0.090), Old Jack Pine (0.120,0.095), Old Black Spruce (0.135, 0.070). Moreover, a mid-arctic summer atmospheric model and a continental aerosol model were selected to characterized the atmosphere above the BOREAS sites. These assumptions may induce slight errors in the procedure, which are not expected however to modify deeply the magnitude and shape of the resulting BRDF.

C.  Results

The plots of either the principal or the perpendicular plane shown in this section were modelled discretely at the same view and solar angles as the measurements. This method of representation was chosen since the aircraft was not flying exactly in the perpendicular or principal solar plane. Continuous BRDF simulations can be seen in Figs.14 and 15.

For the comparison between the model and the measurements, the regression factor (coefficient of determination) R² and the root mean squared (R.M.S.) difference are used as quantificative guide. R² is generally an indication of the curves similarity and the R.M.S. indicate the numerical accuracy. Table 2 has the R² and R.M.S. values for all the simulations. The convention used here is that the backscattering side is plotted with negative view zenith angles. For the perpendicular plane, the half plane that was closer to the backscattering side has been assigned to the negative q_v.

1) Old black spruce

Using the parameters of the tower site, the spatialized dataset is compared with the model in Fig. 9 . Each grid point (every 5 degrees) in Fig. 9(a) and Fig. 9(b) is a weighted interpolation of the measurements. The interpolation averages the measurements and removes some effects of the inhomogeneity of the forests. Each grid point corresponding to a specific q_v(0^° to 60^°) and f (0^° to 360^°) is then simulated using an averaged q_s of 33.5^° to give Fig. 9(c) and Fig. 9(d). For both bands, the model and the measurements have the same general shape except in the hotspot region. The aircraft flight did not allow measurements to be taken close enough to the principal plane so that the maximum in reflectance was not reached. There is a double peak in the measurements because in theory, it is symmetrical on the two sides of the principal plane. Using this symmetrical assumption the data were duplicated on the opposite site, giving a mirror effect. Figures 9e and 9f compare the POLDER measurements (x-axis) and the model (y-axis). Each modelled reflectance was done using the q_v, q_sand f from the measurements. The diagonal lines represents the one-to-one relationship between the model and the measurements. For the red band, a regression factor R² of 0.97 and a R.M.S. difference of 0.0022 are computed. As in Fig. 8a, the lower reflectance values are the ones that are not well simulated, the model overestimating the reflectance on the forwardscattering side. The near-infrared simulation is also very good with a R² of 0.99 and a R.M.S. of 0.0058. The lower reflectance values are also overestimated by the model. The regression factors are close to unity indicating that the shape of the simulations are in good agreement with the measurements.

2)  Young jack pine

The model mimics the observation closely in both red [Fig. 10(a) and10(b)] and NIR [Fig. 10(c) and 10(d)] bands. In Fig. 10a the model did not reproduce the measured hotspot, but it must be noted that the largest measured value is usually not exactly at the hotspot because none of the discrete q_v is equal to q_s. The red band principal plane gives a regression factor of R² = 0.97 with a R.M.S. of 0.002. For the perpendicular plane, R² is 0.85 and the R.M.S. is 0.004. The NIR simulation is generally in good agreement with a R² equal to 0.98 and a R.M.S. of 0.01 for the principal plane. A R² value of less than 0.01 was found for the perpendicular plane but the R.M.S. is also 0.01. This means that the curve's shape is not well modelled, in the perpendicular plane for the NIR band. The two curves diverge on the backscattering side, but the values are very similar. Based on measurements that showed that most shoots were between 10 and 30 degrees [4], a value of 20 degrees (a_L = 80^°) for the shoots main axis from the zenith was used in the simulation. The branch inclination was set at 15^° from horizontal (a_b=15^°). It was observed that the YJP simulations were only slightly improved by using the branch architecture, e.g., without the branch architecture (using G(q) = 0.5 and W_E = 0.72), the regression factor was 0.90 for the red band principal plane.

Figure 11 has the plots for the spatialized dataset compared with the model. The shape of the red reflectance distribution is well modelled but since the difference, although small (R.M.S. is 0.0094), is not systematic, it gives a small regression factor R² of only 0.57. The measurements are more scattered than the simulation points. A comparison of the spatialized data near the principal plane in Fig. 11  and the POLDER reflectance measurements in the red band at the BOREAS tower site (Fig. 10) shows that the latter are in the lower values of the spatialized data . The reason may be that the extent of the forest around the site is limited (400-500 m), so some pixels included in the analysis may be of other cover types (mostly clear-cuts with higher red reflectivity). The NIR simulations also closely reproduce the measurements with a R² of 0.92 and a R.M.S. of 0.012. NIR reflectances for the jack pine forest and the surrounding may be similar.

3)  Old jack pine

2.3.4  Old Aspen

It is shown in Fig. 13a and 13b that the BRDF shape of both planes in the red band BRDF are well reproduced, with R² of 0.93 and 0.82 for the double-canopy and 0.94 and 0.82 for the single-canopy. The R.M.S. diferences between the double-canopy simulation and the data are very good with 0.0027 for the principal plane and 0.0013 for the perpendicular plane. Sligthly larger R.M.S. were found for the single-canopy version with 0.0029 and 0.0039 for the principal plane and perpendicular plane, respectively.

In the NIR band both versions of the model can reproduce very well most view angles except for large q_v values on the backscattering side (Fig. 13 c). An underestimation at view angles larger than the sun angle is found. This phenomenon can not be simulated by the 4-Scale model without accurate consideration of the multiple scattering effect. The regression factor R² are 0.89 and 0.88 for the double and single versions, respectively. The perpendicular simulation gives R² of only 0.37 and 0.31 but the R.M.S. difference are about the same as for the principal plane, i.e., 0.01 to 0.02.

Overall, both versions performed well, and the single-canopy version gives larger reflectance values than the double-canopy version. This is probably due to the absence of the second canopy; the single-canopy induces less shadow, thus giving a larger reflectance.

IV  Discussion

The YJP site was well modelled with or without the branch architecture; the resulting curve were similar in both cases. This is in agreement with the fact that even with shoots oriented mainly between 10 and 30 degrees from the zenith, the YJP does not exhibit an erectophile behaviour as much as in the OJP []. The low foliage density of the OJP requires a better representation of the crown architecture to better model the light and view penetration. The larger LAI in the YJP prevents the architectural structure to influence the overall penetration.

A very large W_s was used for the YJP site to avoid very sharp hotspots, suggesting that at this site, sub-branch architecture (groups of shoots) may be more important in the estimation of mutual shadows within tree crowns than the individual shoots. Large W_s induces less clumped canopy based on gap size theories [8] which comfirms that the clumping of foliage in branches is less important for that site.

The addition of the branch architecture improved the accuracy of the old jack pine site modelling. Figure 14 shows how the branch architecture changes the BRDF of the jack pine site in the near-infrared band. The angle distribution of the shoot (20^° from the zenith) induces more foliage visible at large zenith angle, and with a higher reflectivity for the foliage than for the ground, it increases the NIR reflectance. The shoots orientation also decreases the canopy gap fraction near nadir, but at those angles the canopy view path is done through multiple layer within the same tree crown. Near nadir, the branch architecture does not change much the proportion of background and canopy viewed (see the gap fraction simuation of Fig. 3). One effect of the branch architecture is to clump the foliage into branches; the difference in the simulated BRDF curves between the case with and without the architecture is small because the branch clumping effect has already been partly considered in the use of W_E) for the whole stand.

The effect of the Neyman grouping index (m₂) on the reflectance was not much investigated in the paper. Measurements from the old jack pine site revealed a m₂ of 3, but the sampling area was too small to really see the grouping effect at large scales. The use of more quadrats can show the importance of the Neyman grouping, but it also has the disadvantage of having smaller quadrats which can cause computation errors in the gap fraction at large zenith angles when the view line or the solar beam penetration goes beyond the quadrat from the top of the canopy to the ground surface. The Neyman grouping has a strong effect on the gap fraction, but this effect is in part compensated when the total reflectance is computed after the reflectivities is multiplied by the different scene components. Reflectance from forest with a large contract between the background and canopy reflectivity would be more affected by more clumping of trees.

The mutual shadowing of leaves within tree crowns is modelled with Q_tot1 and Q_tot2. Although the results are good, it seems that they tend to make the hotspot too sharp, especially at large q_s (see Fig. 13(c), and also Fig. 14 in [10]). Figure 15 shows simulations of the principal plane BRDF of a black spruce forest. Without the mutual shadows between leaves (which is equivalent of having Q_tot1 = 1 and Q_tot2 = 0) we see that the overall shape would be much different. At large view angles, on the backscattering side (negative q_v), the modelled BRDF does not show a decreasing trend after the hotspot and is too large for q_v between -40^° and 0^°. This is often seen at large q_s because of multiple-scattering. In general red band BRDF is better simulated then NIR band BRDF, especially near the hotspot. This indicates that the 4-Scale model can be further improved by considering multiple-scattering effects. The model assumes trees of equal height and therefore underestimates the amount of shaded crowns seen. This underestimation could conterbalance the overestimation of shaded foliage computed with Q_tot1 and Q_tot2. 4-Scale needs to be further developed to consider the enhanced mutual shadows in stands with large tree height variation.

In the OA site, the number of hazelnut stems was not an important factor in the resulting BRDF especially when we put enough hazelnut silhouette to cover the ground (3200 stems/hectare and more). The aspen trees are definitely the dominant factors in the BRDF signature of this site. This explains why the BRDF was also well modelled with the single-canopy version of the model. The probability of seeing the canopy illuminated P_T, or P_T1 in the case of the double-canopy version, has the general shape of the BRDF curve. The aspen site provides a challenge because of the complexity of the double-canopy feature, but the model is able to simulate it with a reasonable accuracy. Althought both single and double-canopy version of the model gave similar results, the single-canopy needed to consider the hazelnut reflectivity as part of the background. The OA site also has more multiple-scattering than the other sites, and the angular independence of the multiple-scattering effect may cause considerable errors in this case. For such a complex canopy, a more accurate mathematical description of multiple-scattering is needed.

V.  Conclusion

The POLDER instrument has proven to be a very efficient instrument capable of capturing the signatures of different forest canopies. It has the angular and spatial resolutions required to study the hotspot effect. The quality of the measurements is consistent, although some uncertainties still exist in atmospheric correction at large view angles.

With its multiple input parameters, the 4-Scale model is a flexible tool to study how canopy architecture at each scale influences the reflectance. The reason for the good simulation using only a small parameters set lies in the underlying physics of the interaction between radiation and the canopies described in the model. The parameters and the canopy architecture considerations allows sufficient flexibility of the model to simulate various canopy types. We stress, however, the advantage of this complex model over simpler geometric-optical models in its ability to simulate the BRDF over the entire angle range even with fixed general parameters.